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Abstract 

Max-stable processes have recently been used in modeling spatial extremes. Due to 
unavailable full likelihood, inferences have to rely on composite likelihood constructed 
from bivariate or trivariate marginal densities. Since max-stable processes specify the 
dependence of block maxima, the data used for inference are block maxima from all 
sites. The point process approach or the peaks over threshold approach for univariate 
extreme value analysis, which uses more data and is preferred by practitioners, does 
not adapt easily to the spatial setting. To use more data in spatial extremes modeling, 
we propose a composite likelihood approach that fuses block maxima data with daily 
records. The procedure separates the estimation of marginal parameters and depen- 
dence parameters in two steps. The first step estimates marginal parameters using an 
independent composite likelihood from the point process approach with daily records. 
Given the marginal parameter estimates, the second step estimates the dependence pa- 
rameters by a pairwise likelihood with block maxima. Based on the theory of inference 
functions, the asymptotic consistency and asymptotic normality can be established, 
and the limiting covariance matrix is estimated with a sandwich variance estimator. In 
a simulation study, the two-step estimator was found to be more efficient, and in some 
scenarios much more, than the existing composite likelihood approach based on block 
maxima data only. The practical utility is illustrated through the precipitation data 
at 20 sites in California over 55 years, and the two-step estimator gives much tighter 
confidence regions for risk measures of jointly defined events. 

Key words and phrases: estimating functions, extreme value analysis, spatial depen- 
dence. 
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Weather and climate extremes such as extreme temperature and precipitation have been 



associated with human health in both non- infectious and infectious diseases (Patz et al. 



2003, 2005). Extremes data are often spatial in nature as data are recorded at a network 
of monitoring stations over time. Large scale atmospheric circulation such as El Nino may 
increase the occurrence of extreme weather events over large regions (e.g., Zhang et al. , 2010 ). 



Rare events that occur at multiple locations within a very short time interval can cause more 
damage, consume more resources, and demand stronger emergency management. To mitigate 
losses and be better prepared, a holistic understanding of the dependence of extreme events 
in a spatial context is necessary. Although univariate extreme value modeling has been 
well developed (e.g., Coles, 2001), spatial extremes modeling has gained sharpened focus 



only recently (e.g., de Haan and Pereira, 2006 



Buishand et al. 


2008 


Padoan et al. 


2010 


reviews are 


Davison et al. 


(2012 


) and 


Bacro 



and Gaetan (2012), with the latter one focusing on spatial max-stable processes. 

Max-stable processes are an important tool for spatial extremes modeling. A max- 
stable process is an infinite-dimensional extension of multivaraite extreme value distribution 
such that its every finite dimensional marginal distribution satisfies the max-stable property 
(de Haan, 1984). Specific parametric models have been applied in practice. Two most widely 



used models are the Smith model, also known as the Gaussian extreme value process ( Smith 



1990), and the Schlather model, also known as the extremal Gaussian process (Schlather 



2002). Because the joint density is unavailable when four or more sites are involved, recent 



applications (Padoan et al. , 2010; Davison and Gholamrezaee , 2012) relied on the compos 



ite likelihood approach (Lindsay, 1988) based on annual maxima data for inferences. The 



efficiency of the composite likelihood estimator can be improved if triplewise likelihood is 



used in place of pairwise likelihood, as shown by Genton et al. (2011 ) with the Smith model. 
Nevertheless, the trivariate marginal distribution has not been available for other max-stable 
process models. 

For the univariate case, extreme value analysis can take advantage of more information 
than fitting the generalized extreme value (GEV) distribution to block maxima data (e.g., 
Coles, 2001). When daily records are available, two well known approaches are the peaks 



over threshold (POT) approach (Balkema and de Haan, 1974 Pickands 



generally, the point process approach (Pickands, 1971 Leadbetter et al. 



1975) and, more 
1983). They are 



favored by practitioners in hydrology and climate research since they use more information 
in the data and perform better than the block maxima approach 



(e.g. 



Katz et al., 2002 



Tanaka and Takara, 2002). A natural question is, can these approaches be applied in spatial 
extremes modeling to take fuller advantage of daily record? Unfortunately, a max-stable 
process model only specifies the dependence structure of the block maxima and says nothing 
about other exceedances over the thresholds. In fact, as pointed out by Falk and Michel 
(2009), application of the POT approach in multivariate extreme value modeling poses two 
main problems: 1) which distributions describe the exceedances, and 2) how exceedances are 
defined in a multivaraite setting. Investigation of these problems are still lively continuing 
(e.g., Rootzen and Tajvidi 



2006 Falk and Guillou 2008; Falk et al., 2010). 



1 Without resorting to a spatial version of the POT approach, can we make better use of the 

2 data if daily observations are available and improve the efficiency of parameter estimators? 

3 We propose a two-step composite likelihood approach for inferences in max-stable processes 
modeling with spatial extremes that fuses block maxima data with daily records from each 

5 site. The first step estimates marginal parameters using an independent composite likelihood 
for point processes with daily records data. Given the marginal parameter estimates, the 
7 second step estimates dependence parameters by a pairwise composite likelihood with block 
a maxima data. The two-step approach has been studied recently for multivariate models 



4 



6 



9 to overcome the computational difficulty in maximum likelihood estimation (Zhao and Joe 

10 2005 Joe, 2005). A difference between our two-step approach and theirs is that we use 

11 different data in the two steps; the first step uses daily records while the second step uses 

12 annual maxima. With the separate estimation of marginal parameters and dependence 

13 parameters in two steps, the approach is easily grasped by practitioners. As expected, in 
w our simulation study, the two-step estimator shows substantial gain in efficiency compared 

15 to the existing composite likelihood approach that uses only the block maxima data. Such 

16 gains can lead to much narrower confidence intervals for risk measures such as joint return 

17 levels. 

is The rest of the article is organized as follows. The spatial max-stable process model 

19 defined by all univariate marginal distributions and a spatial dependence structure is intro- 

20 duced in Section [2j In Section [3j we present details of the two-step composite likelihood 

21 approach, the asymptotic properties of the estimator, and how to estimate the limiting vari- 

22 ance. A simulation study is reported in Section [4] to assess the efficiency gain of the two-step 

23 approach in comparison to the current practice based on block maxima data only. The pro- 

24 posed method is applied to the precipitation data from 20 sites in California over 55 years 

25 in Section [5j providing more compact confidence regions for joint return levels. Section [6] 

26 concludes with some discussion. 



2 Spatial Max-stable Process Model 



Max-stable processes are extensions of the multivariate extreme value distribution to the 
infinite dimensional case. Without loss of generality, max-stable processes are often presented 
with unit Frechet margins with distribution function F(z) = exp(— 1/z), z > 0. Such max- 
stable processes are known as simple max-stable processes. A simple max-stable process 
Z(s) on a spatial domain D with unit Frechet margins has all finite dimensional marginal 
distributions satisfying the max-stablility property: 



Pr{Z(xi) < kz u Z(x p ) < kz p } k = Pi{Z(xx) < z u 



[Xr 



< z P }, 



for any p sites {x\ 



,} C D, all Z\ > 0, . . . , z p > and k G N. The max-stablity 



property is equivalently specifying that all p-order marginal distributions are multivariate 
extreme value distributions. 

Two classes of spectral representations of max stable processes have been proposed by, 



38 respectively, de Haan (1984) and Schlather (2002), both of which can be put into a unified 
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form (Davison et al. 2012). Let II be a Poisson process on R + with intensity du/u 2 . Let 
W^(x), x G D, £ > 0, be independent copies of a non-negative stationary process W(x) with 
E{W(x)} = 1 for all x e D. Then, 

Z(x) = max£W^(x), x E D, 

is a stationary simple max-stable process. 

Parametric models for max-stable processes can be constructed by choosing different 
random process W(x). For example, a class of rainfall storm model is obtained by letting 
W(x) = g(x — £), where £ is a random point in D, and g is a density function (de Haan, 1984 
Smith[|l~990| ). This model has a nice interpretation that W(x) is the rainfall at location x from 



a storm with center £ and shape g. If g is a normal density with mean zero and covariance 
matrix E, the resulting process Z(s) is the Smith model, also known as the Gaussian extreme 
value process. This model receives a great interest due to its mathematical tractability and 
nice interpretation (|Smith 1990 Coles, 1993 de Haan and Pereira 



2006 Padoan et al. 



2010 Genton et al. , 2011). The bivariate marginal distribution function at two sites x% and 



x 2 can be shown to be 



Pr[Z(xi) < z x , Z(x 2 ) < z 2 ] = exp 



1 / a 1, z 2 
Z\ \ I a z\ 



1 / a 1, z x 

_$ + -log- 

z 2 \2 a z 2 



where $ is the standard normal cumulative distribution function, a 2 = Ax T S _1 Ax, and 
Ax = x\ — x 2 . The density function can be obtained by differentiating the distribution 
function. 

A max-stable process model for spatial extremes is decomposed into two parts: marginal 
distributions and the spatial dependence structure. Let M(s) be the extremal variable at site 
s in a domain D C R 2 . Let X(s) be the covariate vector available at site s. The marginal 
models are GEV distributions: 



M(s) ~ GEV(/x(s) 



a(s 



22 where /x(s), u(s), and £(s) are the location, scale, and shape parameters of the GEV dis- 

23 tribution at site s, respectively. The covariate vector is incorporated into the parameters 

24 through 



X T (s)f3 a , 



25 where g M , g a , and g% are known link functions for /i, a, and £, respectively, and /3 T = 

26 (Pj , (3j , (3j) is the vector containing all marginal parameters. The spatial dependence is 

27 characterized by a max-stable process (MSP) with dependence parameter 9: 

F- l {G s (M(s);P)}~MSP(9), (2) 



4 



1 where F is the distribution function of unit Frechet, and G s (-;(3) is the distribution function 

2 of GEV(/i(s), cr(s), £(s)) with parameter vector (3. If the Smith model is used, for instance, 9 

3 would contain all three parameters in S, (an, on, a^i)- The whole model is characterized by 

4 r] T = (/3 T , 9 T ), which concatenates the marginal parameters and the dependence parameters. 

5 Since the full likelihood of Z(s) is unavailable, the current practice in parameter estima- 
e tion is the composite likelihood approach constructed from bivariate marginal densities (e.g., 



Padoan et al. 2010). The daily records from each site s do not contribute to the estima- 
s tion except the maximum from block (year), which may otherwise increase the efficiency in 

9 marginal parameter estimation and in turn improve the efficiency in dependence parameter 

10 estimation as well. This motivates our two-step approach. 



3 Two- Step Approach 
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Suppose that we observe the full record of each block with block size m at S sites over n 
blocks (e.g. years, months). Let Y Sjt be the whole observation vector at site s in block 
t, s = 1,...,S, t = l,...,n. Let Y a> t,k De the kth observation within block t at site s, 
k = 1, . . . ,m and M Sji be the block maximum. The block maxima data M = {M s j : s = 
1, . . . , S; t = 1, . . . , n} are observations from the model defined in (|lp and §2§. Let X(s) be 
the covariate vector at site s; examples are longitude, latitude, and elevation. For simplicity, 
we assume that there is no temporal dependence; temporal dependence can be handled with 
declustering in practice. 

Our two-step composite likelihood procedure estimates marginal parameters in the first 
step and the dependence parameters in the second step. 



22 Step 1 The first step is based on an independent likelihood constructed from the point 

23 process approach for univariate extreme value analysis, utilizing all full block data but ig- 

24 noring the spatial dependence. Let Y = {Y s t : s = 1, . . . , S; t = 1, . . . , n}. Let u s be the 

25 threshold chosen for site s, s — 1, . . . , S. The independent loglikelihood has the form 



h(P;Y) 



n S 

EE 

t=l s=l 



If. 



S (AY, 



s,th 



(3) 



26 where 



-Its 



(# y. 



S,t I 



+ 



l+is 

E 

k:Y at ,k>u a 



Us 



O-.s 



- log a s 



-V6 



i Hog i+e 



a s 



t=l l U,s: 



is simply log- 



27 The contribution to the independent loglikelihood from site s, 

28 likelihood of the point process approach in a univariate extreme value analysis (Smith 
1989). Since we assume indepdence from block to block, the contribution from block t 

J2s=i^u,s- The maximizer of ([3]), j3 n , is the estimator of (3. 



30 IS l\t 



5 



1 Step 2 The second step uses only block maxima and estimates dependence parameters 

2 based on a pairwise composite likelihood, with marginal parameters replaced by the estimates 

3 from the first step. The dependence parameters are estimated by maximizing a pairwise 

4 likelihood based on block maxima with /3 fixed at /3 n . Let fij(-', 9, /3) be the bivariate marginal 

5 density of the MSP in ^ at site i and j, with dependence parameters 9 and marginal 
e parameters (3. Define a log pairwise composite likelihood 

n 

l 2 {6- j3 n , M) = ^t{0] L M a>t :s = l,...,S), (4) 
t=i 

7 where the contribution from block t is 

5-i S 

£ 2t (9; 0,M s>t :s = l,...,S) = J2Yl lo S M( M M> M ^ ^ ■ 

i=l j=i+l 

8 Our estimator for 9, 9 n , is the maximizer of Q. 

9 The asymptotic properties of the two-step estimator r)J = (/3j,0j) can be derived 



10 

11 



with general theory of estimating functions (Godambe, 1991). Let ipu0) = diu/df3. Let 



i>2t{P, 9) = d£ 2t /d9. The estimator f\ n is the solution to the estimating equations 

n 

£>fo) = 0, (5) 



t=i 

12 where 

Mv) = 



MP) 
MM), 

13 Under mild regularity conditions, as n — > oo, the solution fj n to ^ is consistent to the true 
w parameter vector i] Q , and 

Vn(77„-7fo) ^N(0,Q), 

15 where Q = A~ l B(A~ 1 ) T is the inverse of the Godambe information matrix, with A = 

16 linv+oo n~ l J2t=i diptiv) / dr] T , and B = lim^^ n' 1 J2t=i ?Pt(v)^7 (v) ■ 

17 With independent replicates at the block level, S can be easily estimated with the sample 
is versions of A and B. We estimate A with 

A If / d^ lt 0n)/d(3 T \ 
n n \dML n )/d(3 T dML 9 n )/d9 T J ' 

19 where the triangular form comes from the fact that the first step does not involve the depen- 

20 dence parameters. This matrix involves the second-order derivatives of the log composite 

21 likelihood, which can be difficult to obtain in practice. An alternative formula that invloves 

22 only the first-order derivatives of the log composite likelihoods is given in the Appendix. For 

23 B, we estimate it with 



B 



1 n 



MMl0n) MMl0nA) 
n \ML On)lpZ0n) ML O n )ML 9 



6 



Our estimator of Q is then Cl n = A^B^A, 



-l^T 



The estimation of B in this problem is 



2 straightforward because we have replicates, unlike in an usual spatial data setting where there 



^s only on e replicate and some bootstrap procedure would need to be developed (Heagerty 



and Lele, 1998 Heagerty and Lumley, 2000). 



4 Simulation Study 



6 



A simulation study was conducted to investigate the performance of the two-step approach 
7 using daily record in comparison with the existing pairwise composite likelihood approach 
s using block maxima only. It also assessed the validity of the sandwich variance estimator of 

9 the asymptotic variance of the estimators from the two-step approach. The study region is 

10 confined to be [— 5,5] 2 . Factors of our experiment are: the marginal model, spatial depen- 

11 dence structure, the number of sites S, and the sample size n (usually corresponding to the 

12 number of years in a real dataset). 

13 The marginal distribution of block maxima at each site s is a GEV distribution with 
location /z s , scale a s , and shape £ s . Let Xi(s) and X 2 (s) denote the latitude and longitude 



14 



16 



15 of site s. We considered first a model of the form 

Us = Py.fi + P^iX^s) + (3y :2 X 2 (s), 

is — /%o> 

where /3 Mi o = 5, p^i = —0.5, (3^2 = 1, Ar,o = 2.5, and /%o = 0.2. The Smith model was 

l? used for the spatial dependence structure. Two dependence levels were considered: weak 

is anisotropic dependence with E = Ei = 4Q and moderate anisotropic dependence with 

19 E = E 2 = 1QQ, where Q is a bivariate correlation matrix with correlation coefficient 0.5. To 

20 make the simulation realistic, we considered the number of sites S G {25, 49} and the sample 

21 size n 6 {20,50, 100}. For each S, the sites form a square grid over [— 5,5] 2 (Figure [T]). 

22 For each scenario, 1000 datasets of daily observations were generated. Note that the daily 

23 observations need to be generated with care such that the annual maxima follow the max- 

24 stable process models specified by the marginal GEV models and the Smith dependence 

25 structure. For each day, we generated a realization from the simple Smith model with 

26 specified dependence level and divided the realization by 365. Assuming independence from 

27 day to day (in a real data analysis declustering would be done to account for temporal 

28 dependence), the componentwise maximum vector at all sites over 365 days follow the simple 

29 Smith model by the max-stability property. The daily series at each site were transformed 

30 to have the specified GEV marginal distributions for the annual maxima. 



31 For a given dataset, the R package Spatial Extremes (Ribatet, 2011) was used to obtain 

32 the pairwise likelihood estimator based on annual maxima data. These estimates were used 

33 as the starting values in the two-step approach, where the general purpose optimizer optim 

34 in R was used to maximize the composite likelihoods in the two steps. The threshold in the 

35 first step of two-step approach was chosen to be the 95th sample percentile at each site. 
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Figure 1: Sites generated for simulation study. Left panel: S = 25; right panel: S = 49. 



1 The empirical mean squared error (MSE) of the estimators from the two-step approach 

2 and pairwise likelihood approach based on the 1000 replicates are summarized in Table [Tj 

3 As expected, the MSEs from both approaches decrease as sample size n increases for both 

4 marginal and dependence parameters. As the dependence level gets stronger from Ei to 

5 S 2 , the MSEs from both approaches increase, and the magnitude of the change is more 
e for dependence parameters than for marginal parameters. The number of sites has little 
7 affect on the MSEs for both approaches. For a given combination of (n, S, £), the two-step 
s approach always gives smaller MSE than the pairwise likelihood approach. The difference is 

9 more substantial for marginal parameters than for dependence parameters, especially with 

10 the regression coefficients in the location parameter, in which cases, the MSE of the two-step 

11 estimator is about 2% of that of the pairwise likelihood estimator. 

12 For a better inspection, Table [2] reports the relative efficiency of the two approaches for 

13 each parameter — the ratio of the empirical MSE from the two-step approach (Mi) to that 
u from the block maxima pairwise likelihood approach (M 2 ). Overall, the efficiency gains of 

15 the two-step approach are large, with the relative efficiency varying from 1.7% to 87%. The 

16 two-step approach improves substantially for the marginal parameters. In particular, we 

17 observe a striking efficiency gain (with relative efficiency between 1.7% and 3.9%) in /^i 
is and /3 M; 2, the coefficients of latitude and longitude in the location parameter. It is also of 

19 interest to look at /^o, since the shape parameter £ governs the tail behaviour of the GEV 

20 distribution and plays an important role in predicting return levels. The relative efficiency 

21 for is between 16% and 24% for all scenarios, which also indicates substantial gain. For 

22 the dependence parameters, the efficiency gain is appreciable with relative efficiency ranging 
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Table 1: Empirical MSE for model 1 based on two approaches {M\\ the existing pairwise 
likelihood approach using block maxima data; M 2 : the two-step approach using data fusion). 



n = 20 n = 50 n = 100 



Si Y*2 Si Y*2 Si S2 



s 




25 


49 


25 


49 


25 


49 


25 


49 


25 


49 


25 


49 


011 


Mi 


1.30 


1.09 


38.94 


44.02 


0.52 


0.45 


16.56 


16.92 


0.24 


0.22 


7.88 


8.39 




M 2 


1.00 


0.82 


29.60 


34.54 


0.39 


0.32 


12.42 


12.91 


0.18 


0.16 


5.89 


6.01 


012 


Mi 


0.70 


0.60 


22.48 


21.92 


0.28 


0.22 


8.05 


8.60 


0.12 


0.11 


3.87 


3.86 




M 2 


0.60 


0.51 


19.04 


18.67 


0.24 


0.18 


6.87 


7.33 


0.11 


0.10 


3.27 


3.08 


022 


Mi 


1.34 


1.18 


50.07 


46.51 


0.54 


0.42 


16.85 


16.52 


0.22 


0.22 


7.84 


7.98 




M 2 


0.97 


0.86 


40.06 


35.96 


0.40 


0.30 


12.94 


13.16 


0.17 


0.16 


5.88 


5.65 


/W xl0 ~ 2 ) 


Mi 


7.49 


7.85 


16.27 


16.73 


2.68 


3.23 


7.47 


7.38 


1.43 


1.42 


3.28 


3.37 


M 2 


5.72 


5.85 


11.27 


12.55 


2.24 


2.34 


5.13 


4.99 


1.05 


1.12 


2.28 


2.28 


&,i( xl0 ~ 4 ) 


Mi 


18.11 


19.49 


26.94 


29.68 


7.10 


7.61 


11.78 


12.00 


3.41 


3.92 


4.81 


5.31 




M 2 


0.40 


0.43 


0.91 


0.54 


0.16 


0.18 


0.22 


0.22 


0.08 


0.09 


0.10 


0.12 


/?M,2(X10- 4 ) 


Mi 


17.74 


20.31 


28.56 


31.11 


7.00 


7.11 


10.46 


11.53 


3.55 


3.73 


4.71 


5.58 




M 2 


0.39 


0.42 


1.11 


0.53 


0.16 


0.16 


0.20 


0.23 


0.08 


0.08 


0.10 


0.10 


/3.,o(xlO- 2 ) 


Mi 


4.02 


4.29 


9.70 


9.95 


1.65 


1.74 


4.04 


3.91 


0.79 


0.80 


1.86 


2.08 




M 2 


3.20 


3.36 


6.84 


7.22 


1.33 


1.29 


3.05 


2.83 


0.59 


0.63 


1.38 


1.40 


/3 C ,o(xlO- 3 ) 


Mi 


4.84 


4.64 


8.89 


8.58 


2.03 


1.88 


3.15 


3.15 


0.96 


0.92 


1.64 


1.72 




M 2 


0.78 


0.83 


1.84 


1.81 


0.35 


0.33 


0.72 


0.68 


0.18 


0.19 


0.39 


0.39 
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Table 2: Relative efficiency in MSE for the pairwise likelihood approach using block maxima 
data relative to the two-step approach using data fusion under marginal model 1. 



71 
1 b 


V 




U 11 


Cl2 


<J 22 


B n 


B i 


B o 


B n 


Be n 


20 


Si 


25 


0.77 


0. 


85 


0.73 


0.76 


0.022 


0.022 


0.79 


0.16 






49 


0.75 


0. 


85 


0.73 


0.75 


0.022 


0.021 


0.78 


0.18 




E 2 


25 


0.76 


0. 


85 


0.80 


0.69 


0.034 


0.039 


0.71 


0.21 






49 


0.78 


0. 


85 


0.77 


0.75 


0.018 


0.017 


0.73 


0.21 


50 


Si 


25 


0.75 


0. 


84 


0.73 


0.84 


0.022 


0.024 


0.80 


0.17 






49 


0.72 


0. 


82 


0.71 


0.72 


0.023 


0.023 


0.74 


0.17 




s 2 


25 


0.75 


0. 


85 


0.77 


0.69 


0.018 


0.019 


0.75 


0.23 






49 


0.76 


0. 


85 


0.80 


0.68 


0.018 


0.020 


0.72 


0.22 


100 


Si 


25 


0.76 


0. 


87 


0.75 


0.74 


0.024 


0.022 


0.76 


0.19 






49 


0.75 


0. 


85 


0.74 


0.79 


0.023 


0.022 


0.78 


0.20 




s 2 


25 


0.75 


0. 


84 


0.75 


0.70 


0.021 


0.022 


0.74 


0.24 






49 


0.72 


0. 


80 


0.71 


0.67 


0.022 


0.018 


0.67 


0.23 



from 71% to 85%. The gain here can be explained by the fact that the marginal parameters 
are estimated more precisely with the daily data in the first step. 

To examine the bias and the performance of the sandwich variance estimator of the two- 
step approach, we summarize in Table [3] the average of biases, the empirical standard error, 
and the average of the standard error from the sandwich variance estimator based on 1000 
replicates. To save space, only results for n e {20, 50} are reported; the case of n = 100 is 
omitted because the results are already good at n = 50. The biases are very small compared 
to the true values for all parameters. The average standard errors are slightly smaller than 
the empirical standard errors for all parameter estimates with n = 20. Consequently, the 
empirical coverage rates of 95% confidence intervals are slightly smaller than the nominal 
rate for some parameters. Nevertheless, with sample size n = 50, the agreement between 
the empirical standard errors and the average standard errors improves, and the empirical 
coverage rates of 95% confidence intervals are reasonably close to 95%. 

We also considered a second marginal model of the form 

Us = #u,o + P^iXi(s) + (3^ 2 X 2 (s) 
a* = &,o + Ax,i x i( s ) + P*,2X 2 {s) 

where ^ = 5, 0^ = -0.5, (5^ 2 = 1, = 2.5, 0^ = 0.2, B^ 2 = -0.2, and /% = 0.2. 
This model had scale parameter depends on the covariates, and generated data with much 
higher variation than the first model. The relative efficiencies in MSE are summarized in 
Table |1J ranging from 22% to 100%. The results indicate that the two-step approach is 
superior to the pairwise likelihood approach for all parameters in this model too, especially 
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Table 3: Bias, empirical standard errors (ESE), average standard errors (ASE) from sandwich 
estimators, and empirical coverage prabability (CP) of 95% confidence intervals for the model 
parameters (Par) based on 1000 replicates with marginal model 1. 

n = 20 n = 50 

£ S Par bias ESE ASE CP bias ESE ASE CP 



Si 25 



S 2 25 



Cll 


0. 


.08681 


0. 


.994 


0. 


.852 


12 


0. 


.07298 


0. 


.768 


0. 


.673 


C"22 


0. 


.08737 


0. 


.984 


0. 


.854 


P«o 


0. 


.00128 


0. 


.239 


0. 


.231 




0. 


00033 


0. 


.006 


0. 


007 




0. 


00016 


0. 


006 


0. 


.007 




-0. 


.00527 


0. 


.179 


0. 


.171 




-0. 


,00821 


0. 


.027 


0. 


.025 


o"n 


0. 


.05451 


0. 


.904 


0. 


.803 


°" 12 


0. 


06876 


0. 


.712 


0. 


629 


0"22 


0. 


09252 


0. 


.922 


0. 


.819 


Pu 


-0. 


00901 


0. 


.242 


0. 


.241 


il 


0. 


00001 


0. 


.007 


0. 


.007 


ft, 2 


0. 


00005 


0. 


.007 


0. 


007 


/^tr 


-0. 


.01379 


0. 


.183 


0. 


.177 


&0 


-0. 


.00964 


0. 


.027 


0. 


026 


0"11 


0. 


.81744 


5. 


382 


5. 


.072 


°" 12 


0. 


63898 


4. 


319 


3. 


.862 


0"22 


1. 


.31173 


6. 


.195 


5. 


.229 


P//,0 


-0. 


01700 


0. 


335 


0. 


.347 


^,1 


0. 


00018 


0. 


010 


0. 


008 


0*2 


0. 


00029 


0. 


Oil 


0. 


007 


&,0 


-0. 


.02125 


0. 


.261 


0. 


.259 


/%0 


-0. 


.01214 


0. 


.041 


0. 


,039 


o-ii 


0. 


.92824 


5. 


806 


5. 


,116 


^12 


0. 


.45541 


4. 


299 


3. 


,755 


0"22 


0. 


92646 


5. 


.927 


5. 


.105 




-0. 


,00363 


0. 


.354 


0. 


361 




0. 


00002 


0. 


.007 


0. 


008 


/3,,2 


0. 


,00031 


0. 


.007 


0. 


008 


Po-,0 


-0. 


,01033 


0. 


,269 


0. 


269 


/%o 


-0. 


.01041 


0. 


.041 


0. 


040 



89. 


.1 


0. 


.01726 


0, 


626 


0. 


,562 


91. 


.1 


91. 


A 


0. 


03012 


0, 


.488 





,441 


92. 


.5 


89. 


.4 


0. 


03365 


0, 


630 


0. 


,566 


91. 


.0 


93. 


,5 


-0. 


.00407 


0. 


.150 





,147 


93. 


.3 


95. 


.9 


0. 


00006 


0. 


.004 


0. 


.004 


96. 


1 


95. 


.1 


-0. 


00005 


0. 


004 





.004 


95. 


.3 


92. 


.5 


-0. 


00432 


0. 


.115 





.109 


92. 


.1 


89. 


9 


-0. 


00624 


0. 


.018 





.016 


91. 


.3 


89. 


.2 


-0. 


00909 


0. 


.568 





.536 


92. 


.5 


89. 


3 


-0. 


.00275 


0. 


.425 





.411 


93. 


.3 


90. 


.2 


-0. 


.00873 


0. 


.544 


0. 


.539 


93. 


.7 


93 


,5 


0. 


00409 


0. 


.153 





.153 


94. 


.8 


95. 


.0 


0. 


00019 


0. 


004 





,004 


95. 


.1 


94. 


.6 


0. 


00003 


0. 


004 





,004 


96. 


.5 


91. 


,7 


-0. 


00284 


0. 


.114 


0. 


.114 


93. 


.7 


89. 


.9 


-0. 


00667 


0. 


.017 





,017 


91. 


.4 


91. 


.8 


0. 


.61117 


3. 


.473 


3. 


.353 


94. 


.5 


91. 


.8 


0. 


.25649 


2. 


609 


2. 


.478 


94. 


.0 


91. 


,4 


0. 


51189 


3. 


562 


3. 


.338 


93. 


,2 


94. 


,7 


0. 


.01510 


0. 


.226 


0. 


.221 


94. 


.8 


94. 


.8 


0. 


.00014 


0. 


.005 


0. 


.005 


94. 


.8 


95. 


6 


-0. 


00034 


0. 


.004 





.005 


95. 


.6 


91. 


,9 


0. 


.01115 


0. 


.174 


0. 


,167 


94. 


.5 


90. 


,5 


-0. 


,00467 


0. 


,026 


0. 


,025 


93. 


.5 


91. 


.5 


0. 


50904 


3. 


,558 


3 


,311 


93. 


.8 


90. 


,8 


0. 


,29361 


2. 


693 


2. 


.428 


92. 


.3 


90. 


.1 


0. 


,54582 


3. 


,588 


3 


.329 


93. 


.1 


93. 


.5 


0. 


00360 


0. 


.224 





.226 


95. 


.3 


95. 


,7 


0. 


00003 


0. 


005 





.005 


94. 


.7 


95. 


3 


0. 


00000 


0. 


005 





.005 


94. 


.8 


92. 


,4 


0. 


00066 


0. 


.168 





.170 


94. 


.7 


91. 


.2 


-0. 


,00627 


0. 


.025 





,026 


93. 


.7 
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Table 4: Relative efficiency of MSE of the pairwise likelihood approach using block maxima 
data relative to the two-step approach using data fusion under marginal model 2. 



71 
1 b 


V 




U 11 


C~i2 


0"22 


B n 


B 1 


B o 


B n 


B 1 


B o 


Ben 


20 


Si 


25 


0.74 


0.73 


0.71 


0.91 


0.63 


0.65 


1.00 


0.35 


0.33 


0.23 






49 


0.76 


0.73 


0.68 


0.85 


0.58 


0.63 


0.91 


0.31 


0.31 


0.23 




E 2 


25 


0.69 


0.74 


0.79 


0.82 


0.66 


0.71 


0.88 


0.44 


0.45 


0.26 






49 


0.75 


0.74 


0.76 


0.84 


0.65 


0.67 


0.86 


0.42 


0.43 


0.27 


50 


Si 


25 


0.75 


0.82 


0.77 


0.96 


0.62 


0.63 


0.95 


0.36 


0.32 


0.22 






49 


0.73 


0.73 


0.68 


0.90 


0.63 


0.67 


0.92 


0.34 


0.35 


0.23 




s 2 


25 


0.75 


0.82 


0.79 


0.87 


0.69 


0.74 


0.94 


0.49 


0.50 


0.29 






49 


0.75 


0.76 


0.73 


0.77 


0.60 


0.60 


0.84 


0.43 


0.41 


0.26 


100 


Si 


25 


0.81 


0.87 


0.76 


0.87 


0.57 


0.55 


0.93 


0.34 


0.30 


0.25 






49 


0.75 


0.84 


0.76 


0.88 


0.61 


0.65 


0.89 


0.31 


0.32 


0.24 




s 2 


25 


0.77 


0.83 


0.78 


0.80 


0.68 


0.65 


0.85 


0.47 


0.41 


0.29 






49 


0.74 


0.79 


0.72 


0.82 


0.70 


0.68 


0.83 


0.49 


0.42 


0.28 



for the shape parameter (20-29%) and the regression coefficients of latitude and longitude 
in the scale (31-50%) parameter. In contrast to model 1, the efficiency gain for /3 M) i and B^ 
of model 2 is not as large as that of model 1, This may be due to the simplicity of model 1, 
in which latitude and longitude only appear in fi s . 



5 Extreme Winter Precipitation in California 

Daily precipitation records at all monitoring stations in California were extracted from 
the second version of the Global Historical Climatology Network (GHCN), compiled and 
quality controlled at the National Climatic Data Center of the National Oceanic and At- 
mospheric Administration (available online at http://www.ncdc.noaa.gov/oa/climate/ghcn- 



daily/). The same data source was used in Shang et al. (2011). The raw data contains daily 



rainfall observations at 230 sites in California starting from 1878. As most precipitations in 
California occur in winters, only the winter extreme precipitation is considered. The winter 



season here spans from December 1st to March 31st in the following year (Zhang et al. 



2010). Due to missing data, the maximum observation in a given year at a given site was 
considered as the annual maximum only if there were no more than 5% missing daily obser- 
vations during the winter season. In our analysis, we used a balanced data set consisting of 
daily winter precipitation from 1948 to 2002 for 20 sites (see Figure [2j. 

We propose a spatial max-stable process model for the California data. Three explanatory 
variables are available for the marginal parameters. Let X T (s) = {Xi(s) , X 2 (s) , X 3 (s)} be 
the longitude, latitude, and elevation at site s, respectively. The first two are in degrees, and 
the elevation is in 100 meters. They are centered at the values at San Francisco (—122.38, 
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Figure 2: Locations of the 20 monitoring stations in California. The three sites in solid 
circles are Napa, Winters, and Davis. 



37.62, 0.02). The marginal GEV model at site s has parameters 



= 0m,o + ^,i x i( s ) + P»,2X 2 (s) + f3^ 3 X 3 (s) 
cr(s) = p ff>0 + ^ili(s) + /3 a>2 X 2 (s) + P a , 3 X 3 (s) 



This is the same model as that in Shang et al. (2011) except that the Southern Oscillation 
Index term is removed, which otherwise would introduce temporal nonstationarity at all 
sites. The spatial dependence structure is modeled by a Smith model with dispersion matrix 
E. 

To estimate the parameters with the two-step approach, we chose the 95th sample per- 
centile as the threshold u(s) for each site s. As in standard univariate extreme value analyses, 
there was temporal dependence in daily precipitation and we needed to remove clustering 
from the observed data before applying the two-step approach. A simple way of declustering 
is to define consecutive exceedances of a threshold to belong to the same cluster; the cluster 
is terminated once an observation falls below u(s) (e.g., Coles, 2001). 



Table [5] summarizes the parameter estimates and their standard errors from the two- 
step approach and the pairwise likelihood approach. Both approaches gave qualitatively the 



same results, which are consistent with those in Shang et al. (2011). As the covariates were 



centered at the values at San Francisco, the intercepts for the location and scale parameters 
are interpreted as the location and scale, respectively, for the GEV distribution in San 
Francisco. The estimated shape parameter are positive and significantly away from zero, 
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Table 5: Summaries of parameter estimates and the standard errors for two approaches (M^ 
composite pairwise likelihood approach using block maxima data; M 2 : two-step approach 
using data fusion). 



Mi 



Mo 



Model Parameter 


Estimate 


Standard Error 


Estimate 


Standard Error 


Marginal GEV 










Location 


Intercept 


4.677 


0.167 


5.692 


0.073 




Latitude 


— U. lyo 


U.U4y 


— U. lol 


U.Uzo 




Longitude 


— U.o t o 


U.Uoo 


— u. ( y i 


n noo 
U.Uoz 




Elevation 


U.loU 


U.U1 ( 


U.loo 


U.U1 ( 


Scale 


Intercept 


2.307 


0.127 


2.372 


0.069 




Latitude 


-0.145 


0.052 


-0.129 


0.024 




Longitude 


-0.358 


0.067 


-0.358 


0.031 




Elevation 


0.061 


0.018 


0.059 


0.011 


Shape 


Intercept 


0.173 


0.030 


0.092 


0.016 


Dependence structure 












on 


0.766 


0.085 


0.301 


0.036 




012 


-1.382 


0.153 


-0.494 


0.061 




0"22 


2.538 


0.274 


0.882 


0.106 



indicating that the tails of the marginal GEV distributions are heavier than that of a Gumbel 
distribution. Sites in the west and south, with higher elevation tend to have higher location 
parameter and higher scale parameter (higher variation). The spatial dependence among 
extreme precipitations is captured by the parameters in S. The estimated a±i is about a 
third of the estimated a 2 2, indicating anisotropy in the strengths of the spatial dependence. 
The 0"i2 is significantly negative, suggesting a rotation in the ellipsoidal dependence along 
the coast line. 

The two-step approach yielded marginal parameter estimates that are of similar size 
to those from the pairwise likelihood approach except for the shape parameter; the shape 
parameter estimates are 0.173 and 0.092 from the pairwise likelihood approach and the 
two-step approach, respectively. The dependence parameter estimates from the two-step 
approach are about a third to a half of those from the pairwise likelihood approach. For 
this application, We observe more difference in dependence parameter estimates than in 
marginal parameter estimates between the two approaches. As far as the standard errors 
are concerned, the two-step approach gave much smaller standard errors than the pairwise 
likelihood approach except for one parameter — the coefficient of elevation in the location 
model. Most of the reductions are about a half or more. 

What are the implications of the reduced standard error in inferences? Of course, reduced 
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Table 6: Joint 50-year return levels (cm) for three pairs based on both approaches (Mi: 
pairwise likelihood approach using block maxima data; M 2 : two-step approach using data 
fusion). 



Pair 


M 1 




M 2 




95% CI 


Width 


95% CI 


Width 


Napa & Winters 


(8.84, 10.82) 


1.98 


(10.63, 12.17) 


1.54 


Napa & Davis 


(8.83, 10.78) 


1.95 


(10.62, 12.20) 


1.58 


Winters & Davis 


(10.84, 13.61) 


2.77 


(12.34, 14.03) 


1.69 



standard errors in marginal parameters lead to more precise inference about marginal risk 
measures such as return level for each site. Since our context is spatial extremes modeling, 
it is of interest to see how the reduced standard errors affects risk measures of jointly defined 
events. The first one is the joint 50-year return level for two sites s\ and 82, which is defined as 
y 50 such that Pr (Y(si) > y^, Y(s 2 ) > y^o) — 1/50. It is exceeded simultaneously at the two 
sites once every 50 year. Because the bivariate marginal densities of the max-stable process 
are known, y^ can be found numerically for any given parameter vector. The second one 
is the joint sampling distribution of the maximum of the extremal precipitation over every 
50 years at multiple sites. From the max-stability property, the dependence structure of the 
joint distribution is characterized by the same max-stable process up to a scale. Realizations 
from the distribution can be drawn for any number of sites and can be used to assess the 
joint behavior. 

Consider three stations near the Sacramento area: Napa (122.25°W, 38.27°N), Winters 
(121.97°W, 38.52°N), and Davis (121.78°W, 38.53°N). We generated model parameters from 
the multivariate normal approximation of the estimator, with a large number (N = 5000) of 
realizations. For each generated parameter vector, the joint 50-year return level was obtained 
numerically for each pair of sites, since the bivariate marginal distribution function at two 
sites is known for the Smith model. Table [6] shows the 95% confidence intervals of joint 50- 
year return levels for three pairs with estimator from both approaches. As can be seen, the 
two-step approach predicts higher joint return levels than the pairwise likelihood approach 
but the intervals still overlap. The two-step approach always gives a tighter confidence 
interval, about 61-81% of the lengths from the pairwise likelihood approach. In particular, 
the joint return level for Winters and Davis are much higher than those for the other two 
pairs, which could be explained by the closeness in distance and thus strong dependence 
between them. The two-step approach gives an interval 39% percent narrower that the 
pairwise likelihood approach for this pair. 

For all the three sites, we now look at the joint sampling distribution of the maximum 
of every 50-year period. For each of the 5000 parameter vector drawn from the asymptotic 
normal distribution, we generate one realization for the 50- year sitewise maxima. Figure [3] 
shows the scatter plot of the 5000 draws of 50-year maxima for each pair based on both 
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Precipitation(cm| Precipitation(cm) Precipilation(cm| 



Figure 3: The 50-year sample return levels (cm) for three pairs on the log scale. Upper: 
pairwise likelihood approach using block maxima data; Lower: two-step approach using 
data fusion. Left: Napa & Winters; Center: Napa & Davis; Right: Winters & Davis. 

1 approaches on the log scale. The two-step approach always gives a more concentrated scatter 

2 plot than the pairwise likelihood approach, thus giving a more accurate prediction. There 

3 seems to be a positive dependence between each pair, which is especially obvious for the last 

4 pair (Winters and Davis), because of their strong dependence. 

5 6 Discussion 

e In contrast to the existing composite likelihood approach which utilizes only block maxima, 

? our two-step approach fuses block maxima data with daily records and makes more efficient 

s inferences about the parameters. The first step estimates marginal parameters based on 

9 an independent marginal likelihood from the point processes approach for univariate ex- 

10 treme value analysis; the second step estimates dependence parameters using a pairwise 

11 composite likelihood by block maximum data, with marginal parameters replaced by their 

12 estimates from the first step. This method is simple, intuitive for practitioners who are 
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familiar with univariate extreme value modeling, avoiding defining multivariate thresholds 
and multivariate Pareto process modeling (Aulbach and Falk, 2012). Our simulation study 



showed substantial efficiency gain from the existing composite likelihood approach to the 
two-step approach with realistic number of sites and sample sizes. In an application to 
extreme winter precipitation in California at 20 sites over 55 years, the two-step estimator 
gives much smaller standard errors for all parameters as well as tighter confidence intervals 
for joint risk measures than the pairwise likelihood approach. 

Further research is merited in several aspects to extend the method. In the first step, we 
did not address threshold selection, an important and still active problem even for univariate 



extreme value analysis (e.g., Guillou and Hall, 2001 Thompson et al. , 2009). Recent research 
showed promising approach with quantile regression for nonstationarity with covariate in- 



formation (Northrop and Jonathan, 2011). Alternatively, we may use the r largest order 
statistic to construct the marginal likelihood (e.g., Coles, 2001). This way, there is no need 
to specify the threshold and temporal nonstationarity may be incorporated to detect change 
in trend. In the second step, for max-stable process models that have known trivariate 
marginal densities, a triplewise likelihood can be used in place of the pairwise likelihood; for 



the Smith model, efficiency gain has been shown by Genton et al. (2011). More broadly, an 



iterative approach may be developed by adding a likelihood term involving the dependence 
parameters to the marginal likelihood in the first step to borrow strength from other sites 
using the spatial dependence in estimating the marginal parameters. How to appropriately 
weight the pieces in the composite likelihood and its practical utility need to be investigated. 
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24 A Alternative Formula for A 



n 



25 We first look into the contribution of each site and each pair at each year to the terms in 

26 A n . Let = de lt Jdp. Let = ^Sfij{(M i>u M j>t ); 9, 0) and ^,^(0,9) = 

27 d£2t,(i,j) 1 99 . Then, A n can be rewritten as 



A, 



1 n 

n |J V E?j=iH<j fy*,W)(Pn, 9 n )/dp H, d^^, 9 n )/d6 



28 Instead of calculating the second-order derivatives in A n , we use the first-order derivatives 

29 based on the second Bartlett identity, assuming that the univariate and bivariate marginal 

30 models are correctly specified. Let (p2t,(i,j)W, 9) = <9^2t,(i,j)/<9/3. We can estimate A by 
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